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Abstract 

The numerical simulation of colliding solitary waves with compact support 
arising from the Rosenau-Hyman K(n, n) equation requires the addition of 
artificial dissipation for stability in the majority of methods. The price to 
pay is the appearance of trailing tails, amplitude damping, and delays as 
the solution evolves. These undesirable effects can be corrected by prop- 
erly counterbalancing two sources of artificial dissipation; this procedure is 
designed by using the slow time evolution of the parameters of the solitary 
waves under the presence of the dissipation determined by means of adia- 
batic perturbation methods. The validity of the tail removal methodology 
is demonstrated on a Pade numerical scheme. The tails are completely re- 
moved leaving only a small compact ripple at the original position of their 
front, and the numerical stability of the scheme under compacton collisions 
is preserved, as shown by extensive numerical experiments for several values 
of n. 
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1. Introduction 



Rosenau and Hyman [11] came across solitary waves with compact support, 
therein referred to as compactons, while studying the nonlinear dispersion in 
the formation of patterns in liquid drops. These traveling waves arose as 
solutions of a particular generalization of the well-known Korteweg-de Vries 
equation, commonly denoted as K(m,n), that reads 



u t + 



{u m ) x + 



0. 



m > 0, 1 < n < 3, 



(1) 



where u(x,t) is the wave amplitude as a function of the spatial variable x, 
and time t. Current analytical studies of this equation 0, 0, 0, 0, 0, 0, 0, 0], 
cannot deal with the interaction between compactons. 
methods, such as pseudospectral schemes fl L 10[, finite elements |lll.ll2l 



differences lll|, |13| . Pade approximants 



Hence, numerical 
finite 
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or 



modified equations 

particle schemes based on the dispersive- velocity method [13], must be used. 

The majority of the numerical methods for compacton equations require 
the addition of artificial viscosity in order to cope with collisions, since oth- 
erwise instabilities appear which may cause a blow up in the solution. In 
pseudospectral schemes a hyperviscosity term is used where a second-order, 
linear dissipative term affects only the high-frequency modes of the solu- 



tion by using a high-pass filter (lOj. In both finite element and finite dif- 



ference schemes a fourth-order, linear dissipative term is frequently incorpo- 



rated [111 . Il2l . Il3l Il5l . Il6| . The addition of artificial viscosity in Eq. ([T]) distorts 
the original compacton, generating small trailing tails, amplitude damping 
and velocity losses. The development of a method for the numerical stabi- 
lization of the compactons in collisions without the appearance of tails is an 
open problem to be addressed in this paper. 

Numerically-induced phenomena can be studied and corrected by means 
of the method of modified equations 16j, [l8(. For the analysis of these phe- 



nomena, perturbation methods [19f| can be applied to the analysis of effects 
introduced by local truncation errors as perturbations of the original evo- 
lution equation. These methods have been successfully applied to ordi nary 



differential equations [18|, |20j and nonlinear evolution equations [2ll . l22l |23 



For compacton equations, adiabatic perturbation methods have been con- 
sidered for the K(n, n) equation in Garralon and Villatoro |24| . extending 
previous results for the n = 2 case in Pikovsky and Rosenau [25|, and Rus 
and Villator o pisf - These techniques, which generalize previous results for 
solitons 27|, |28|, |29j, can be used to analytically determine the trailing tails 
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introduced by the artificial viscosity, opening the possibility of their removal 
by means of correction terms modifying the original evolution equation. 

In this paper adiabatic perturbation methods are applied to the K(n, n) 
equation and a new numerical technique for trailing tail removal is intro- 
duced. The validity of this procedure is checked by means of a Pade nu- 
merical method. Section [2] is devoted to the development of the tail removal 
technique based on the adiabatic perturbation to the numerical scheme with 
artificial viscosity. In section [3] we briefly present the numerical method based 
on Pade approximants used to solve the K(n, n) equation. The results of ex- 
tensive numerical experiments are presented in Section H] for several values 
of n. Final conclusions are summarized in Section [SJ 

2. Adiabatic perturbations 

The adiabatic perturbation method is applied to the K(n, n) equation 
with both second- and fourth-order linear dissipations. Such a method de- 
termines the slow time evolution of the parameters of the compactly sup- 
ported solitary waves by using that of the invariants under the dissipative 
perturbation. 

Let us consider the perturbed K(n, n) equation given by 



where 1 < n < 3, the perturbation V(u) is a function of u and its spatial and 
temporal derivatives, and \e\ <C 1 is a small parameter. After multiplying 
Eq. (j2J) by u n and integrating in space, the only non-null term in its left-hand 
side is the first one, resulting in 



whose left-hand side is the temporal derivative of the second invariant of 
the unperturbed K(n,n), i.e., exactly nil for Eq. (EJ with e = 0. However, 
perturbations such that the right-hand side of Eq. ([3]) is non-zero result in 
the adiabatic evolution of parameters of the compacton solution of the un- 
perturbed equation. Introducing the slow time scale r = et, the compacton 
solution of the K(n, n) can be written as 




(2) 





2nc(r) 
n + 1 




n 



1 



) 



l/(n-l) 



u c (x,t,r) 



2n 



(x-c(r)t) 



(4) 
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for \x — c{r)t\ < n7c/(n — 1), and u c (x,t,r) = otherwise. Inserting this 
ansatz into the perturbed equation ([3]) the slow time evolution of the velocity 
c(t) of the perturbed compacton can be calculated; note that the amplitude 
of the perturbed compacton is uniquely determined by this velocity. 
Let us determine the evolution of c(r) for perturbation given by 

eV{u) = a 2 u xx - a 4 u xxxx , (5) 

where |ct2 1 , |cf4 1 *C 1 are small parameters; this perturbation is dissipative if 
a 2 > and a 4 > 0. The substitution of Eqs. (EJ) and fll]) into Eq. yields 
an ordinary differential equation for c(r) written as 

(n-1) 2 (n -l) 3 ((n -3) 71-1) 

c (t) = -, tt«2 c r -. — 5- — a 4 c(r). (6) 

w n (n + 3) w (n - 5) n 3 (n + 3) w w 

The solution of this equation is 

< x / m / (n - 1 ) 2 (n-l) 3 ((n-3)n-l) \ 

c M = c(0) exp (-^j a, r - ( „ ! ^,,3 ( ,/ + 3) «. rj , (7) 

showing that, under the perturbation, the velocity, as well as the amplitude, 
decays and produces lags compared with the unperturbed compacton. 

The perturbation ([5]) introduces trailing tails in the perturbed compacton 
not accounted for by the ansatz (0J, since the perturbed equation does not 
possess solitary waves solutions; hence the amplitude of the initial pulse 
slowly decays as described by Eq. ([7]). This decay results in the formation of a 
plateau behind the compacton; this behaviour is similar to that of the solitons 



of the Korteweg-de Vries equation under the same kind of perturbation [28 



Let us note that the adiabatic perturbation method can be used to estimate 



the shape of this tail, as shown in Refs. [25] and [26] for the K(2 } 2) equation, 
and Ref. |24] for the K(n, n) equation. 

The main contribution of this paper is the introduction of a new proce- 
dure for removing the trailing tail by properly adjusting the values of the 
small parameters a 2 and 0:4, suggested by the fact that their contributions 
in Eq. ()5]) have opposite signs. Fixing the value of a 2 as a function of a 4 and 
n yields 

/ x (n — 1) ((n — 3) n — 1) 

For 1 < n < 3, the value of 02(11, a^)ja^ is positive, increasing from 1/9 to 
1/4 as n decreases from 3 to 2, and decreasing from 1/4 to as n decreases 
from 2 to 1. 
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3. The numerical method 



The tail removal procedure presented in the previous section could be in- 
corporated into a numerical method for Eq. ([1]) which uses artificial viscosity 
in order to deal with compactons collisions. Obviously, this procedure could 
introduce new instabilities. For illustration purposes, let us consider one of 
the most used methods for Eq. (pQ), a Pade approximation in space, with 
periodic boundary conditions in a finite interval [0, L], and a method of lines 



in time as described in [14J, [15|. Equation (jSJ) with perturbations (jSJ) and a 



moving frame of reference with velocity Cq is given by 

Eq[u) = - c u x + (u n ) x + (u n ) xxx - a 2 u xx + a 4 u xxxx = 0. (9) 
The Pade method of this equation results in 

A(E)^> - coB(E)^ + B(E)(C/ i ) n + C(E)(U,) n 

-a 2 S(E)U j + a i V(E)U j = 0, j = 0,l,...M, (10) 

where the spatial grid nodes are Xj = j Ax, for j — 0, 1, . . . M, with Ax = 
L/M, the numerical solution is Uj(t) ~ u(xj,t), the shift operator E is de- 
fined as EUj = U j+1 , and A~ 1 (E)B(E), ^- 1 (E)5(E) ^- 1 (E)C(E), and 
.4. _1 (E) £>(E) are Pade operators for the first-, second-, third-, and fourth- 
order spatial derivatives respectively. Concretely, the Pade operators to be 
used in this paper are 

E" 2 +26E" 1 +66 + 26E 1 + E 2 
A(E) = 



and 



B(E) = 
5(E), 
C(E) 

2>(E) = 



120 

-E^-IOE^+IOE^E 2 
24Ax : 
E~ 2 +2E~ 1 -6 + 2E 1 + E 2 
6Ax 2 ' 
- E~ 2 +2 E" 1 -2 E 1 + E 2 
" 2Ax 3 ' 

E" 2 -4E _1 +6 -4E X + E 2 
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corresponding to the following approximations to the first four spatial deriva- 
tives 

A-\E)B(E)u( Xj ,t) = ^(x 3 ,t) + ^ ^( Xj ,t) + O (Ax 8 ) 
A- 1 (E)S(E)u(x j ,t) = Q^(xj,t) + ^ ^{xj,t) + O (Ax 6 ) 



and 



A-\E)C{E)u{x^t) = g(x„i) - ^ ^(x,,t) + O (Ax 6 ) 



Ax 6 


d 7 u 


5040 dx 7 


Ax 4 


d G u 


720 


dx 6 


Ax 4 


d 7 u 


240 


dx 7 


Ax 2 


d 6 u 


12 


dx 6 



A- 1 (E)V(E)u(x j ,t) = ^(si,f) - + (Ax 4 ) . 

The modified equation of the numerical method ( ITUl) . which includes the local 
truncation error terms, is given by 

Eq[u{ Xj ,t)] -<*< — ^(x„t) + O (Ax 4 ) = 0, (11) 

showing that the numerical scheme is dissipative (for a 4 > 0) after the ap- 
plication of the tail removal procedure. 

For the discretization in time of Eq. (fTUl) . let us use the second-order 
accurate, implicit midpoint rule, yielding 

A{E) 3 - — j -(c B(E) + a 2 S(E)-a 4 V(E)) 3 - - 3 

B(E) - C[E) ) ( 3 - i ) =0. (12) 

where At is the time step and C/j ~ u(xj, z At). The resulting nonlinear sys- 
tem of equations is solved for f/j i+1 ^ by using the Newton's iterative method. 

4. Presentation of results 

Let us study the behaviour of the new procedure for trailing tails removal 
on the compacton solutions of the K(n, n) equation with n G {3, 2, 5/3, 3/2, 
7/5, 4/3, 9/7, 5/4}. Numerical solutions with several values of will be 
compared with both a 2 = and a 2 = a 2 (n,a4), for the propagation of a 
one-compacton solution and for the collisions between two compactons. 
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4-1. One-compacton solutions 

Let us first consider the evolution of a one-compacton solution with and 
without the tail removal procedure. Figure [1] shows snapshots of the numer- 
ical solution at the same instant of time for several values of n. The left 
plots show the trailing tails that appear when oli = 0, for = 10~ 3 (top 
plot) and 0:4 = 10~ 5 (bottom plot). Note that a zoom in has been used since 
the amplitude of the tails is several orders of magnitude smaller than that 
of the compacton. The tail starts at the initial location of the compacton at 
t = with a ripple with a negative peak followed by a positive one; after the 
ripple, the tail has a nearly constant plateau connecting with the left edge of 
the compacton. The amplitude of the tail decreases as n does and as so does 
«4. The right plots show the removal of these tails when using a^w, 04), for 
04 = 1CT 3 (top plot) and = 1CT 5 (bottom plot). The plateau completely 
disappears, but the ripple in the tail front persists with smaller amplitude. 
In long-time integrations, the compacton is far away and clearly separated 
from the ripple located at its initial position, so it behaves as if there is no 
tail. 

Figure [2] compares the evolution in time of the compacton maximum 
amplitude with «2 = (dashed lines) and ac2(n, 0:4) (solid ones), i.e., with 
and without the tails. The appearance of the tails compensates the loss in 
the compacton maximum amplitude due to the dissipation since the first 
invariant is preserved during the propagation, hence, with the tail removal 
procedure, the maximum amplitude practically retains its initial value. 

Table [1] shows the difference between the location of the analytical and 
numerical compacton maximum amplitudes at a time t = 2000 for several 
values of a^, with a 2 = and 0:2(^,04). The numerical velocity of the 
compactons in the dissipativeless case (when a 2 = a 4 = 0) is very near, but 
not exactly equal, to the analytical one, resulting in a very small difference 
in the position at t = 2000 of about 0.6 (for the values of Ax and At used 
in Table [1]), which is mostly independent of the value of n. The use of the 
tail removal procedure, 02(^,04), results in exactly the same value for this 
difference in location at t = 2000, except for either n = 3, or a 4 = 10 -2 
and n > 2, as shown in Table [TJ Without the tail removal procedure, the 
difference between the numerical and analytical compactons after 20000 time 
steps can be very large, in fact, larger than 26, 3, and 0.8 (in units of space) 
for 0:4 = 10~ 2 , 10~ 3 , and 10 -4 respectively. Hence, Table [I] shows that the 
tail removal procedure works properly and results in a numerical velocity 
equal to that of the numerical method without dissipation. 
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Figure 1: Numerical one-compacton solution of the K(n,n) equation for n = 2, 5/3, 3/2, 
7/5, 4/3, 9/7, and 5/4, with c = 1, Co = 0.5, Ax = 0.05, and At = 0.04. The top plots use 
an = 10~ 3 and the bottom ones 0:4 = 10~ 5 ; the left plots use a 2 — and the right ones 
0(2{n, 014). 
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Figure 2: The evolution in time of the maximum amplitude of the one-compacton solution 
of the K(n,ri) equation for n = 2, 5/3, 3/2, 7/5, 4/3, 9/7, and 5/4, with c = 1, a 4 = 1CT 3 , 
c = 1, Ax = 0.05, and At = 0.04. The dashed lines use a 2 — and the solid ones 
a.2(n, a 4 ). 





«4 


= lO" 2 


0:4 


= 10~ 3 


0:4 


= 10" 4 


n 


a 2 = 


a 2 (n, a 4 ) 


a 2 = 


a 2 (n, a 4 ) 


«2 = 


a 2 (n, a 4 ) 


3 


505.50 


108.50 


53.20 


5.50 


5.80 


0.90 


2 


404.90 


0.90 


47.80 


0.60 


5.40 


0.60 


5/3 


223.10 


0.70 


25.20 


0.60 


3.10 


0.60 


3/2 


129.80 


0.60 


14.40 


0.60 


2.00 


0.60 


7/5 


80.80 


0.60 


9.00 


0.60 


1.40 


0.60 


4/3 


53.20 


0.60 


6.00 


0.60 


1.10 


0.60 


9/7 


36.70 


0.60 


4.20 


0.60 


0.90 


0.60 


5/4 


26.30 


0.60 


3.20 


0.60 


0.80 


0.60 



Table 1: Difference between the maximum amplitude location of the analytical and nu- 
merical one-compacton solutions at time t = 2000 for = 10~ 2 , 10~ 3 , and 10 -4 , with 
0L2 = and a2(n, a±). In this table, c = 1.0, cq = 0.5, Ax = 0.1, and At = 0.1 have been 
used. 
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4-2. Compacton collisions 

The numerical simulation of compacton collisions requires the use of arti- 
ficial viscosity in order to avoid the appearance of instabilities. It seems that 
the tail removal procedure introduced in this paper may affect the stability of 
the simulations for colliding compactons since the explicit artificial viscosity 
terms has been canceled. Let us summarize our extensive set of simulations 
showing this is not the case. 

Figure |3] illustrates the collision of two compactons of the If (5/3, 5/3) 
equation with velocities c\ = 1.0, the tallest one, and c 2 = 0.5. The upper 
left plot shows both compactons before the collision. Two snapshots of the 
collision are presented in the upper right and the lower left plots. Finally, 
the lower right plot indicates that both compactons recover their amplitude 
and velocities, but a small ripple remains at the location of the interaction. 
This collision is typical for all the K(n, n) equations. 

Figure H] shows a zoom in of the numerical simulation of the collision 
between two compactons of the K(2, 2) equation. The initial positions of 
the compactons and their velocities (ci = 1, c 2 = 0.5 and c = 0.1) have 
been selected in order to ensure that the fronts of the tail of each compacton 
are far away at the time of collision. Left plots show both compactons for 
«o = 0, before (top plot) and after (middle and bottom ones) their collision. 
The amplitude of the tail of each compacton depends only on the value of 
a 4 (equal to 10~ 3 in this figure) when the tails are well separated (see upper 
left plot in Figure H|) , but when the fastest compacton gets to the trailing tail 
of the slowest, the tail of this second one rides over that of the first one (se e 
middle left plot). As it is widely known in compacton literature (ll.[l2. 13. 30 



a ripple appears after the collision which can be seen around x = 500 in the 
lower left plot; the ripple amplitude is about two orders of magnitude larger 
than that of the tails. 

The application of the tail removal procedure cancels the trailing tails 
without affecting the stability of the numerical method during the collision 
of the compactons, as illustrated in the right plots shown in Fig. H] for the 
K(2, 2) equation. The small ripples located at the front of the tails which 
is not removed by taking a 2 (n, a 4 ) (see upper right plot) do not affect the 
collision which occurs as if they were absent. Lower right plot illustrates that 
the compactons, the ripple associated to the collision, and those associated to 
the tails behave as compactly supported solutions (except for the numerically 
induced, self-similar, backward and forward radiations reported by Rus and 



Villatoro [3J| , whose amplitude is smaller than 10 in the plots of Fig. H]) 



10 
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650 450 
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Figure 3: Four snapshots of the collision of two compactons of the if (5/3, 5/3) equation 
with velocities c\ = 1.0 and C2 = 0.5 at times t = 370, 380, 400, and 450, from upper left 
to lower right. The parameters used in this simulation are ct2 = 0, 0:4 = 10~ 3 , cq = 0.1, 
Ax = 0.1, and At = 0.1. 
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Figure 4: Numerical simulation of two colliding compactons of the K(2, 2) equation with 
velocities c\ = 1.0 and C2 = 0.5 before (top plots with t = 150) and after (middle and 
bottom ones with t = 700) their mutual collision, with ai = (left plots) and 02(^,04) 
(right ones). The parameters used in this simulation are = 10~ 3 , cq = 0.1, Ax = 0.1, 
and At = 0.1. yi 
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Figure 5: Numerical simulation of two colliding compactons of the -ftT(5/3, 5/3) (top plots) 
and K(7/5, 7/5) (bottom ones) equations with velocities c\ — 1.0 and c 2 = 0.5 before (left 
plots) and after (right ones) their mutual collision, with a2(n, oti). The parameters used 
in this simulation are ay = 10~ 3 , c = 0.1, Ax = 0.1, and At = 0.1. 
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The results shown in Figure H] for the K(2, 2) equation are representative 
of an extensive set of simulations of compactons collisions for n G {3, 2, 5/3, 
3/2, 7/5, 4/3, 9/7, 5/4}, different values of a^, ci, c?, Cq, Ax, and At. Figured 
shows representative results for the if (5/3, 5/3) (top plots) and K(7/5, 7/5) 
(bottom ones) equations. The slowest compacton has been stopped by using 
Co and the plots show that the tails have been removed, except for the ripple 
at their fronts, which are two orders of magnitude smaller than the residual 
after the collision of the compactons. Further results for other values of n 
are omitted here for the sake of brevity. In all the cases the proper behavior 
of the tail removal procedure has been observed. 

In long-time integrations, using periodic boundary conditions, two com- 
pactons collide multiple times and their trailing tails ride over one another 
increasing their total amplitude. In such simulations, the value of the artifi- 
cial viscosity must be chosen properly in order to avoid the appearance of 
instabilities resulting in the sudden blow-up of the solution inside the time 
integration interval. Our results show that the tail removal procedure do not 
affect the time of blow-up, validating its good performance in terms of the 
preservation of the stability of the numerical scheme. 



5. Conclusions 

A procedure for trailing tail removal for numerical methods for the K(n, n) 
equation incorporating artificial viscosity has been introduced. This proce- 
dure is based on the analysis of the effect of the artificial viscosity in the 
propagation of compactons by means of the adiabatic perturbation method. 
The performance of the new procedure is illustrated by using a widely used 
numerical method based on Pade approximants (which also can be derived 
from finite element and finite difference formulations). The results obtained 
after an extensive set of simulations show the effectiveness of new the trailing 
tail removal in the propagation of solutions with both only one compacton 
and two compactons in mutual interaction. 

The new procedure could be applied to other numerical schemes for the 
K(n, n) equation and for the numerical study of other nonlinear evolution 
equations presenting compactly supported solutions, such as those with co- 



sine/sine compactons reported in Ref. [32 . 
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